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In thermoacoustic regenerators, the interaction of thermo-viscous boundary layers and 
axial temperature gradients causes a conversion from thermal energy to acoustic power or 
vice versa. In this paper, an improved analytical model for thermoacoustic boundary layer 
effects in the presence of mean flow is derived and analyzed. Previous formulations of the 
thermo-acoustic effect take into account effects of mean flow on acoustic propagation 
only implicitly, i.e. in as much as mean flow influences the mean temperature field. The 
new model, however, includes additional terms in the perturbation equations, which 
describe explicitly the interaction between steady mean flow and acoustics. For a parallel 
plate pore the three-dimensional thermoacoustic equations are derived and reduced to a 
transversally averaged system of differential equations by applying Green's function 
technique and suitable assumptions. The resulting one-dimensional perturbation equa¬ 
tions are then solved numerically for two sets of boundary conditions to obtain the linear 
scattering matrix coefficients. The solutions, generated for a wide range of frequencies, 
can be applied in a low-order “network model" context to study the stability of ther¬ 
moacoustic devices. The impact of mean flow on the thermoacoustic interaction is 
investigated and validated against full computational fluid dynamics simulations of 
laminar, compressible flow for one specific configuration. It is shown that at low 
frequencies (Womersley number < 1) the new formulation predicts the acoustic behavior 
more accurately than the earlier formulations. Finally, the ideas and benefit of further 
improved and more complex models for higher Mach numbers are discussed. 

© 2014 Published by Elsevier Ltd. 


1. Introduction 

The interest in the interaction of acoustics and heat flow dates back more than 200 years [1,2], One of the investigated 
issues is the thermoacoustic (TA) effect: in narrow geometries thermo-viscous boundary layers interact with axial mean 
temperature gradients. This interaction causes a conversion of thermal energy to acoustic power or vice versa. Sondhauss [3] 
observed this effect first. He described the occurrence of spontaneous loud sound emissions in narrow pipes during the glass 
blowing process. Although Rayleigh [4] had already given the first qualitative explanations for this effect, only incomplete 
descriptions of the phenomenon were provided by different authors [5,6], Rott was the first author to completely describe 
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Nomenclature 

0 

blockage ratio 



¥ 

generic variable, axial distribution 

b 

inhomogeneous term 

CO 

angular frequency, rad s _1 

c 

speed of sound, m s~ 1 



c p 

specific heat at constant pressure, J kg -1 K 1 

Subscripts and abbreviations 

Cy 

specific heat at constant volume, J kg -1 K 1 



5 

characteristic wave amplitude, m s _1 

0 

mean quantity 

ts 

characteristic wave amplitude, m s _1 

1 

acoustic quantity 

h 

specific enthalpy, J kg -1 

a 

acoustic 

k 

body force, m 2 s~ 2 

m 

boundary 

I< 

thermal conductivity, W m 1 K 1 

d 

downstream 

L 

length, m 

u 

upstream 

La 

Lautrec number 

lin 

linear 

Ma 

Mach number 

exp 

exponential 

Ma a 

acoustic Mach number 

BC 

boundary condition 

P 

pressure, Pa 

CFD 

computational fluid dynamics 

Pe 

Peclet number 

CFL 

Courant-Friedrichs-Lewy 

Pr 

Prandtl number 

LODI 

local one-dimensional inviscid 

R 

radius, m 

LES 

large eddy simulation 

Re 

Reynolds number 

M/SO 

multiple input single output 

U 

flow velocity, ms~' 

SI 

system identifications 

Wo 

Womersley number 

SISO 

single input single output 

9? 

gas constant, J kg -1 K 1 

GF 

Green's function 

t 

time, s 

ODE 

ordinary differential equation 

T 

temperature, K 

A 

modeled term 

u 

velocity, m s _1 

f 

concerning bulk viscosity 



IMP 

implicit model 

Greek characters 

I< 

thermal 



M I 

model I 

P 

compressibility, I< 1 

M II 

model II 

Y 

heat capacity ratio 

V 

viscous 

S 

boundary layer thickness, m 

ref 

reference 

£s 

influence factor of solid-fluid interaction 

S 

solid 

e 

ratio of dimensions 

TA 

thermoacoustic 

n 

ratio of dimensionless quantities 

E 

energy conservation 

8 

transversal distribution 

G 

gas law 

K 

Helmholtz number 

G 

averaged gas law 

M 

dynamic viscosity, kg m 1 s~ ] 

k 

intermediate averaged continuity 

V 

kinematic viscosity, m 2 s _1 

K 

final averaged continuity 

Q 

density, kg m~ 3 

T 

energy conservation 

o 

thermal conductivities ratio 

X 

axial momentum conservation 

T 

shear stress tensor, N m 3 




the phenomenon analytically in a series of publications [7-12]. Wheatley [13], his successor Swift [14] and finally in't 
panhuis [15] and in't panhuis et al. [16] improved the analytical understanding of the TA effect. Arnott et al. [17] enhanced 
the applicability of this approach for arbitrary pore shapes. Watanabe et al. [18] further improved these models by 
accounting for nonlinear effects such as drag or heat transfer. The TA effect has also been investigated in various 
computational fluid dynamics (CFD) simulations (e.g. [19-23]). Their analytical formulations facilitate the description of the 
axial propagation of acoustic velocity u, and pressure p, in terms of their cross-section averaged quantities. These 
transversally averaged parameters (p,), (Ui) form a set of transport equations in the axial pore direction. Those transport 
equations are applied in a one-dimensional numerical tool [24] to predict the operating conditions in TA devices. 

Previous publications have presented ideas for probable applications of TA devices [13,14], and some authors claim that 
the technology is ripe for the market [25], Nevertheless, almost 30 years after the first attempts of technical realizations of 
such apparatus, the technology is only applied in niche applications [26-28], Thus, the main aim of researchers remains to 
expand the applicability of TA for cooling and power generation. 

All cited authors restrict their descriptions to quiescent mean flow conditions. Until now, only a few publications have 
tried to deal with mean flow affected TA power conversion [29-31], Bammann et al. [32] provide an overview over the 
published literature on flow-through TA refrigeration. One of the main reasons for the modest interest in their topic is the 
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lack of proper modeling of the TA interaction effects under such conditions. All hitherto referenced publications consider the 
impact of mean flow only implicitly by adapting the axial mean temperature profile, or do not address it at all. Based on the 
assumption that the mean flow is due to streaming effects induced by acoustic quantities, a second-order enthalpy flux 
consideration leads to a distortion of the mean temperature profile [8,33-35], The present paper presents and analyzes an 
improved quasi-one-dimensional method, which explicitly takes the laminar mean flow into account in the perturbation 
equations. Using this approach enables the prediction of the performance of TA devices that exploit mean flow as a positive 
effect, e.g. Reid's refrigerator [30], Due to the resulting one-dimensional form, the computational effort for the calculation of 
one configuration is comparatively small. Conceivable applications may, for example, arise in exhaust gas streams, where 
the waste heat enthalpy flux is directly transformed into acoustic power. 

The derivation of the quasi-one-dimensional model is based on the method applied by in't panhuis et al. [16,15], who 
derived TA transport equations that account for almost arbitrary stack pore geometry. Like all other proposed approaches, 
the model takes into account the so-called long-pore assumption, i.e. the hydraulic diameter 2 R is much smaller than the 
characteristic length in axial direction I. Further, in't panhuis [15] follows the common approach of neglecting the 
contribution from mean flow when linearizing the basic set of equations. In the present work, this assumption is not 
invoked and mean flow is explicitly incorporated. 

Taking into account the contribution of mean flow leads to an increase in the number of terms in the considered 
equations. Hence, for the purpose of demonstrating the impact of explicitly accounting for mean flow, the solution is derived 
for a generic problem resulting in the least complex system of transport equations. The geometrical conditions and thermo¬ 
physical correlations are chosen to keep the equations short. Sticking to the general solution with only a few necessary 
assumptions is still feasible; however, this complicates the interpretation of the impact of mean flow due to the huge 
numbers of terms inside the system of equations. Thus, to simplify the model where possible, assumptions will be 
introduced. These assumptions due to compactness of equations are marked as such. 

Apart from spatial distributions of the acoustic quantities for a single frequency, acoustic scattering matrices are used for 
validating the one-dimensional models against CFD results. The latter validation tool is mostly used to investigate the 
stability of the system in the linear regime by using ID network models [36], Although TA engines generally operate in 
nonlinear limit cycles, the identified growth rate at the onset of the acoustic instability is a quantitative measure for the 
achievable pressure amplitudes. Altogether, as the simplicity of the solution is a central request, this publication is restricted 
to the linear acoustic regime. If nonlinearities like the acoustic feedback into the mean temperature distribution or 
saturation effects have to be taken into account, all steps can be processed for the nonlinear situation. 

The paper is organized as follows: In Section 2, the geometrical and mean flow conditions are presented. In Sections 3 
and 4, the one-dimensional acoustic transport equations are derived for a parallel plate stack pore. The flow chart of the 
derivation is sketched in Fig. 1 : First the compressible Navier-Stokes equations (NSE) are non-dimensionalized using narrow 
geometry assumptions. These dimensionless NSEs experience a simplified modal expansion. The mean flow field is solved 
analytically from the resulting zeroth-order set of equations. These solutions are used as input to the first-order set of 
equations, the linearized dimensionless NSEs. This set of equations for the acoustic field is then solved in segregated steps: 
First Green's function (GF) solution technique is applied on the axial momentum and energy equations. This yields 
formulations for their transversal shapes. Closure assumptions have to be made for two convective terms to obtain an 
analytical solution of these formulations. Integrating the resulting equations in the transversal direction creates a set of 
averaged perturbation equations. Finally, eliminating density and temperature perturbations leads to one-dimensional TA 
transport equations in terms of velocity and pressure oscillations. The resulting set of equations accounts for non-zero mean 
flow contributions, which appear explicitly in some terms. Thus, models based on this system of equations are named 
“explicit” in the rest of the paper. The existing model [37] is denoted as “implicit” (IMP), because in that model mean flow 
influences the system of equations only implicitly through its impact on the mean temperature profile, which is computed 
from a second-order energy equation that considers enthalpy transport by the mean flow [8], For the explicit equations two 



Fig. 1. Derivation adapted from in’t panhuis et al. [16], which considers mean flow only implicitly (IMP). In this approach, the mean velocity is explicitly 
taken into account. Model approximations M 1, M II must be invoked to obtain analytical solutions. 
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upstream downstream 

L_ L _J 


Fig. 2. Sketch of the observed problem: The acoustics in a slab pore, i.e. a narrow channel of length L and constant hydraulic radius R within a non- 
isothermal axial field and mean flow are contemplated. The sinusoidal lines indicate the profile of the acoustic quantities p, and u,. These quantities can be 
transformed into characteristic wave amplitudes 5, CO. They form the base for the scattering matrix notation which correlates these incoming and outgoing 
waves at both ends of the channel. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.) 

alternative approaches for the necessary closure assumptions are investigated. Prediction accuracy of these two approaches, 
denoted by M I and M II, is tested for the spatial acoustic propagation for a single frequency and scattering behavior for 
different test cases against a full CFD simulation and the existing implicit one-dimensional solution method. Finally, an 
outlook to improved modeling approaches is presented. 

2. Problem description 

2.1. Geometric settings 

A channel of length L and constant hydraulic radius R is one of the most simple geometries to be considered for a TA 
system. Restricting the problem to a constant pore height simplifies the contact conditions of the fluid and solid regions in 
the problem. Choosing a channel geometry leads to a two-dimensional problem in space, which enormously reduces the 
size of the derived equations. The solid plates in the problem sketched in Fig. 2 are 2 R s thick. The origin of the Cartesian 
coordinate system is positioned at the upstream (u) end in the center of the channel. The counter side of the pore at x=L is 
indicated as the downstream (d) position. Note that any quantity featuring an index S refers to the solid stack. 



2.2. Thermo-physical conditions and mean parameter assumptions 

The solid region inside the pore is assumed to consist of a homogeneous material. All solid properties are only functions 
of the mean temperature. The walls are non-permeable for the working fluid, which is assumed to be an ideal gas. The 
viscosity v, thermal conductivity K and the specific heat capacities c p , c v are assumed to be functions of mean pressure p 0 
and temperature T 0 . The index 0 refers to (temporal) mean quantities, while an index 1 denotes a time dependent 
perturbation quantity. 

The blue to red colorized gradient in the lower part of Fig. 2 indicates the temperature distribution inside the pore. Here, 
the blue indicates the cold side, whereas the hot downstream side is marked in red. 

For the purpose of comparison with experiments, the terms upstream and downstream, which indicate a mean flow 
from the cold to the hot side are introduced. In this paper all reference parameters are chosen to be the cold upstream 
values. The mean flow is presumed to be laminar. Although the incompressibility assumption (see [38]) might be violated by 
the strong impact of the axial temperature, a Poiseuille-like mean flow field is expected. 

3. Non-dimensional, linearized Navier-Stokes equations of a slab pore 

This section covers the first four steps of the sketch depicted in Fig. 1 : At first the Navier-Stokes equations are non- 
dimensionalized. The resulting set of equations is then split into transport equations describing the mean flow field and the 
propagation of acoustic quantities. Here, a special modal expansion technique is applied. 

As boundary layer effects are the cause of the TA effect, the Navier-Stokes equations are a suitable choice of equations 
describing the motion of the fluid under consideration. The equations in the notation of Landau and Lifshitz [39] read as 


D p 

Dt=^ V - V 

(1) 

Dv 

(2) 

= - Vp+V ■ r+pk 

^rg + v . ( Kvn + T: Vv. 

(3) 


where 

Df d'F dT 

-=—=-FV— 

Dt dt dx 
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denotes the substantial derivative of a transported variable *F. The stress tensor t of a Newtonian fluid is described by 

r = 2/4 (vv + (W) T ) +C(V ■ V)I. 

As the body forces k are negligible in the inspected problem, the corresponding term can be omitted. The set of differential 
equations is closed using an equation of state relating density and pressure. Since, in this case, the ideal gas law 

P = P'XT (4) 

is used, the transport of enthalpy h - cf. Eq. (3) - can be rewritten in terms of temperature T 

pc p ^ = /frg+V.(KVD + T:Vv. (5) 


3.1. Non-dimensionalizing the NSEs 

The system of equations (1)—(5), (the NSEs) are non-dimensionalized with respect to the characteristic scales of the 
geometry under consideration. Olson and Swift [40] first derived a set of dimensionless parameters describing TA 
phenomena by applying the Buckingham El theorem. Here, we proceed in the same manner as in't panhuis et al. [16] for 
a two-dimensional set of equations formulated in Cartesian coordinates. The corresponding dimensionless quantities are as 
follows: 

x = Lx, y = Ry, t = — f, u = c Tef u, v = ec ref v 

c ref 

2 

~ 2 ~ -r ^rpf ~ ,, Tp,ref 7, 

d = <?reff?, P=drefCrefP> T=-^-T, C = C ref C, /J = 0 

Cp.ref C ref 

K=K [e fK, c p = C pre fCp, T= ^ f, ^ = // re fp, C = CrefC 


using the index ref for the reference parameters. All axial terms are non-dimensionalized using the length L of the pore, 
while for all transversal terms, v.y and the shear tensor r, the hydraulic radius R of the pore is used. The resulting 
dimensionless numbers are listed in Table 1. The viscous (p), the volume viscous (f) and thermal ( I< for gas, S for solid) 
penetration depths <5are defined as follows: 

2/*ref ^ ref ^ _ I 2/C re f 

(Wq r ef V V 




Table 1 

“/7-Group” for thermoacoustics in narrow pores. 


Symbol 

Formula 

Description 

r 

fp 

Heat capacity ratio 

e 

R 

Ratio of dimensions 


L 


a 

l<s 

Thermal conductivities ratio 


Kref 


<p 

^gas 

Blockage ratio 


A 


K 

coL 

Helmholtz number 

Ma a 

U 

Acoustic Mach number 

Fr 

u 2 

Froude number 


gt 


La 

R 

Fluid Lautrec number 


Sk 


La s 

R 

Stack Lautrec number 




Pr 

s 2 

Prandtl number 


4 


Wo 

R 

Womersley number based on 8 V 




Wo c 

R 

Womersley number based on 8 f 


Si 
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In dimensionless form, the governing equations now read as 

Y~ 1 


P = ' — Tq 
Y 


du dv 
Dt — ^\dx + dy 


Dq 


(7a) 

(7b) 


Du dp 
e Dt = ~- + 


(TU 


dx Wo 2 dy 


-+e 


* 

(tfv tfu \ 

3* d 2 u ' 

Wo 2 

1 dx dy) 

Wo 2 dx dy 


(7c) 


0 = - 


dp 


cPu 


dy wo 2 dy 2 


(7d) 


DT Dp k 
e Di~ /)T Dt + 2La 2 


d 2 T 


dy‘ 


+ 0{e 


(7e) 


dT 


&T 


Qs dt 2La: 


I * 2 


+ 0(e 


(7f) 


From here on, the tilde, characterizing dimensionless variables, is omitted for reasons of clarity and readability. Here, the 
transversal momentum equation (7d) is multiplied by the square of the ratio of dimensions e = R/L Hence, the smallest non¬ 
vanishing term contains information of O(e 0 ). In the equation of state (7a) T, p and q correlate with the specific heat capacity 
ratio y = Cp fc v . The Womersley numbers Wo and Wo- denote the ratio of hydraulic radius to the viscous boundary layer 
thicknesses and <5 f respectively. Their thermal equivalent are the solid (s) and fluid Lautrec numbers La. La and Wo scale in 
terms of the square root of the Prandtl number, i.e. Pr=S 2 /S Although La is generally used to characterize TA devices by 
their porous media [41 J, the authors use the Wo notation because of its familiarity from aero-acoustic applications. The 
frequency dependent component of the Helmholtz number k = q>L/c in the viscous and thermal terms is balanced by the 
square of one of these four dimensionless numbers. The set of transport equations (7f) describes the full fluid dynamics of a 
perfect gas inside a narrow pore at low Ma in a non-dimensional manner. 


3.2. Separating acoustics by modal expansion 

In the next step, the acoustic perturbations of the flow field variables are separated from their mean values. Here, as the 
mean flow velocity is related to the speed of sound c, a special form of modal expansion of the flow variables in terms of the 
acoustic (index a) Mach number Ma a = iq /c re f is applied. This procedure assumes that any variable tqx, t ) may be split into a 
temporally constant mean quantity ToW and a harmonically oscillating part T' 1 (x)e 1Kt . The full standard expansion then 
reads as 


T{\, t) = ‘f'o(x)+Ma a l f / r(x)e ,Kt +0(Ma 2 ). (8) 

In addition, it is assumed that cross-sectional variations in the mean flow quantities are negligible, leading to 

T'o(x)^ f\ 0 (x)dy = ('F 0 )(x). (9) 

Jo 

This implies a simplification, which is exact only for zero mean flow conditions. Its general validity is discussed in the next 
section. 

Olson and Swift [40] showed that an oscillation of the material properties of the size of an acoustic perturbation only 
contributes to higher order terms in the resulting system of equations. Therefore, we consider material properties to be 
functions of the mean quantities, especially leading to 

,J= (To)(x)' (10) 

As we restrict ourselves to linear acoustics, the acoustic Mach number Ma a is considered to be small. The same holds for 
the geometric ratio of dimensions e in the investigated problem. Moreover, we follow the arguments of in’t panhuis [16] 
stating that the ratio of dimensions is of the same order as Ma a . Expressing the latter in terms of q = e 2 /Ma 2 in the set of 
equations (7f) yields individual systems of equations in different orders in e. All expressions have to be satisfied separately 
by the flow. Zeroth-order terms describe the mean flow conditions inside the pore. Terms of first order determine the 
acoustics inside the considered field. As the acoustic Mach number can be chosen arbitrarily, we define i; to be unity. 
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3.3. Mean transport quantities 


The acoustic system of equations contains terms scaling with mean quantities. Thus, finding an analytical description for 
this field implies that the mean parameters can be formulated in terms of functions in x. These formulations are derived 
from the mean transport equations, which are similar to equations (7f) for vanishing e. In addition, the neglect of transversal 
changes in mean flow quantities is discussed. 

A brief look at the transversal momentum Eq. (7d) leads directly to a block profile of the mean pressure and therefore 

PoW = p 0 (x). (11) 

As this publication is restricted to the laminar regime, the mean flow is of Poiseuille type. According to the parallel plate 
solution for this flow type, the transversal velocity v 0 (x) component vanishes 

v 0 (x) = 0. (12) 


Further, all transversal dependency is neglected in order to maintain compactness of the desired equations. Hence, all 
other mean flow quantities are treated to be constant in y, i.e. 'Fo(x) = T'o(x)- It has to be stated that this approach is 
nonessential for T 0 , u 0 and thus q 0 - The consequences of dropping this simplification are discussed in Section 7. This constant 
mean assumption may seem unphysical for the description of the mean flow conditions; however, as we focus on the 
acoustics inside the tube the interactions of the acoustic and mean flow profiles have to be considered. In the implicit 
transport equations (see [14]) some terms depend on the ratio of geometric length to acoustic boundary layer thicknesses. 
The viscous dissipation scales with H„(y) and the thermal relaxation is a function of H«{y) which read as 


cosh((l+i)Woy) _ cosh((l+i)Iay) 

cosh((l+i)Wo) 3n K ~ cosh((l+i)La) 


(13) 


Further, the interaction of both functions, i.e. Hfx) — PrH K (x), affects the TA effect. Considering their profiles with respect to 
the assumed mean flow profile allows a qualitative insight into the interplay of combined terms. The thinner colored lines in 
Fig. 3 give a qualitative impression of their profile for three different Wo = 0.1,1.10 in air. The thicker black lines indicate 
the two parabolic and constant dimensionless mean flow profiles. Due to the friction forces at walls all velocity components 
are small near the wall. For low Wo, the profile of all terms is nearly parabolic, whereas for higher values, the terms show a 
small peak near the wall before they almost reach a constant value far away from the wall. The intersection of the mean 
velocity profiles at y = 2/3 indicates the point from which the contribution of the acoustic profile in interaction with a 
constant mean flow is overestimated with increasing y. 

In the acoustic equations, the product of some mean flow quantities i// 0 (x)</) 0 (x) appears. Considering the product of the 
cross-sectional averages ( y/ o (x))(0 o (x )) as mentioned in the last section instead leads to a certain amount of error. Writing 

V /o(x) = <^o>(x)(l+¥ / o(y)) (14) 


^o(x) = <0 O )(X)(1 + ^oCy)) (15) 

generates four terms for the upper product 

Vo(x)<fi 0 (x) = (iy o )(x)( < p o )(x)(\ + ‘F o (y)+0 o (y)+‘F o (y)0 o (y)). (16) 

If i// 0 denotes the mean flow velocity its transversal contribution is limited by -2/3 < <f'o(y)> < 1/3. In a first approximation, 
the ratio of the hydraulic radius of the pore to the thermal penetration depth 

S T = ---„„ (17) 

C'ref Cp, ref C re fMd 

scales the cross-sectional temperature profile. Ma denotes the reference mean flow Mach number Ma = u 0u /c re{ . In the 
conditions met in TA applications, this ratio is very small. Thus, in contrast to neglecting the second right-hand side term of 
Eq. (16) assuming the other terms to be small induces only small errors. The same argument accounts for terms including q 0 - 
Applying the approach stated in Eq. (15) on the ideal gas law (7a) directly yields the same order of magnitude for the 
contribution of the mean density variation in the transversal direction. Summarizing these two considerations, the error 
committed by the constant mean flow assumption has to be thoroughly investigated and, if necessary, the equations should 
be expanded to accommodate the parabolic profile. 

If the impact of the neglected terms is small, the mean parameter field may be described in terms of cross-sectionally 
averaged quantities. At first, the mean temperature field is considered. If the transversal mean temperature dependency is 
neglected, the conductivity also has to be expressed in terms of a cross-sectionally averaged quantity, which describes the 
combined fluid and solid heat conduction. Inserting the blockage ratio 0 = A gis /A, this mean conductivity (K) becomes 

(K) = 0 -<P)K s + 0K. 

Using this quantity and introducing the Peclet number 

P Cp ref(9 re f C re fL 

Kref 


(18) 

(19) 
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Fig. 3. The three spatial contributions to the acoustic parameters at quiescent conditions according to Eq. (13) for air at Wo = 0.1 (thin solid line), Wo = 1 
(thin dash-dotted line) and Wo = 10 (thin dotted line). A constant (thick dash-dotted line) and a parabolic (thick solid line) mean flow profile are included 
for comparison. For small Wo the spacial contribution resembles a parabolic shape while at larger Wo they show a constant value at some distance to 
the wall. 


for describing the enthalpy flow directly yields 


CpPe(e 0 u 0 ) 


d<T 0 > 


= 7rJ«<) 


dx dx 


A T o) 

dx 


with given boundary values (T)lo = 1 /(y— 1) and (T)|, = T d /[{y-\)T u ]. 
The integral form of the averaged continuity Eq. (7b) reads 


<u 0 <?o>(x) = <u 0 >(0)<eo>(0) = Ma 


( 20 ) 


( 21 ) 


if block profile inflow conditions are met. Inserting Eq. (21) into (20) finally yields 


CpPeMa 


d <T 0 ) 
dx 


d^ 

dx 



( 22 ) 


For isothermal properties, the temperature profile follows an analytic exponential solution. In general, the variation of (K) 
with temperature of real materials is non-negligible and is accounted for in the later computation. These investigations 
showed that the basic nature of the profiles remain quasi-exponential for different ceramics and metals. 

Writing the separate cross-sectional average instead of the left-hand side of Eq. (21) is only exact, if one of the mean 
quantities does not vary in the y-direction. Due to the consideration of transversal averages, the second and third terms of 
Eq. (16) vanish. As discussed above, the cross-sectional variation of the mean density is small compared to its average value. 
Hence, neglecting the product of the y-variations in velocity and density leads to a small underestimation if we assume 

(Q 0 u 0 )(x) = Ma^( eo ){u o>(x). (23) 


Finally, for closure, the axial pressure profile has to be considered. For this purpose, the mean pressure gradient term 
occurring in Eq. (7c) is simplified by considering the Poiseuille flow assumption. Written in a non-dimensional form and 
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assuming a parabolic mean velocity profile leads to 


dp 0 (X) K cf-Up _1_ 

dx — Wo 2 dy 2 Re c e 


(24) 


Further, if the Reynolds number 



(25) 


is larger than the inverse of the ratio of dimensions e, which is often valid for TA applications, the right-hand side term of 
Eq. (24) becomes small. Under those conditions, the influence of the axial mean pressure gradient vanishes. Thus, it is 
assumed that the axial mean pressure profile is constant with respect to its mean quantity. This leads to a pure mean 
temperature dependency of the axial (u 0 ) profile, which is order of magnitudes larger than the changes caused by the mean 
pressure gradient. The impact of the pressure gradient, which physically causes the mean flow, sets the inlet velocity. Using 
this assumption, all mean quantities may be computed. Further, if all the presented assumptions apply, those terms in the 
acoustic system of equations that have contributions in terms of v 0 , the y-derivative of mean quantities, or the axial 
derivative of the mean pressure, may be neglected. 


3.4. Linearized Navier-Stokes equations of a slab pore 


The analytical formulations derived in the last section sufficiently describe the mean flow terms occurring in the acoustic 
set of equations. The acoustic field is governed by the terms of first order in e and in this context they are denoted as the 
linearized Navier-Stokes equations (LNSE) of a slab pore. In the most general form they read as 


di 



(26a) 



d(h do 

(?i+U^- L +T 1 Ul = -<? 
dxdx 


(dv-i diq\ 

V dy dx) 


(26b) 


e 



din 

u 1+ u^ = 

A, 


dp i 
dX 


Kp ifU] 
Wo 2 dy 2 


0 = - 


dfh 

dy 


dTj 


dT 


dT 


iKC p QT i + uqc p — + c p — qUi + c p — uqi = ixrpi + Tu-^-+ 


dp j kK ft 2 !] 


dx 


dx 


dx 2La 2 dy 2 


a 2 


(26c) 


(26d) 

(26e) 




kK s cPTs ,i 
2 Laj dy 2 ’ 


(26f) 


where all first-order quantities are defined in two-dimensional space, i.e. 'F\ = l F](x,y). Flere onward the index 0 as well as 
the spatial averaging <) symbol are dropped for all mean quantities. It is important to note that the gray boxed terms 
explicitly include mean velocity contributions. This is the main difference to previous publications [14,16,29,30]. All these 
terms describe the explicit, spatially averaged impact of mean flow inside a pore. The two-dimensional acoustic quantities 
are affected in- and outside of the thermo-viscous acoustic boundary layers. The system of equations of in't panhuis et al. 
[16] is recovered, if all terms including mean flow effects are dropped. 


4. Transversally averaged acoustic perturbation equations 

The cross-sectional averaging of the LNSEs of a pore (Eqs. (26)) leads to a system of one-dimensional, coupled differential 
equations. This set can then be evaluated numerically with little effort. Nevertheless, the averaged equations contain 
numerous terms. To avoid human errors in this derivation step, the process was carried out with a computational software 
capable of symbolic equation handling [42]. The derivation of the intermediate and final formulations in written form is very 
extensive. The reader is referred to Table 2, containing an overview of each single step of the derivation procedure. The 
integration of the trivial transversal momentum equation (26d) in step 1 yields a constant pressure profile over the pore 
width. In the second step, applying the method of Green’s functions on the axial momentum equation (26c) and modeling 
the term Ai leads to an analytical description of the transversal profile of the acoustic velocity iq. With these results, the 
solution for the energy equations ((26e) and (26f)) is obtained in a similar manner (step 3). The fluctuating temperature in 
the solid domain T 1-S is used to express the non-constant coupling wall temperature T^ in the fluid domain. Inserting the 
expressions of these steps in terms of acoustic velocity Ui, pressure p i and temperature Tj into the gas law (26a) (step 4) 
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Table 2 

Overview of the segregated mathematical steps applied to the LNSEs that lead to a transport ODE system in terms of p, and e. 


Step 

Equations 

Applied methods 

Result 

1 

(26d) 

Transversal integration 

p,(x,y) = p 1 (x) 

2 

(26c) 

(a) Modeling of A^ 

(b) Green's function method 

(c) Transversal integration 

Si=21 1 i(xXuiXx) 

3 

(26e) and (26f) 

(a) Modeling of A 2 

(b) Green's function method 

(c) Boundary coupling 

(d) Transversal integration 

(T 1 )(x)=/«u 1 >(x),<p 1 >(x)) 

4 

(26a) 

(a) Transversal integration 

(b) Inserting 

<diXx)=/«Ui>(x),<Pi>(x)) 

5 

(26b) 

(a) Transversal integration 

(b) Inserting 

(v, >(x) = 0 ^= A 2X (x)(u, >(x) +A u (x)(p l >(x) 

6 


Combination of Steps 2 and 5 

i (;>)■-«(;,) 


yields an expression for the density fluctuations ei- Eliminating all terms except for tq, p lt finally leads to a second ordinary 
differential equation (ODE) of the transversal quantities (step 5) which forms a system of transport equations with the 
solution of step two (step 6). 

A standard analytical method for solving ODEs with a not vanishing source term is Green's functions method. Applying it 
to inhomogeneous differential equations allows one to find an analytical solution for any right-hand side inhomogeneity 
term. The present derivation is based on that technique. Consequently, the general procedure is described first. 

The LNSEs are a set of coupled differential equations containing partial derivatives in all coordinate directions. For the 
purpose of finding an analytical solution in the transversal (y) direction, the relevant acoustic parameters are split into a 
product of ¥ , i(x) = !/q(x)i9(y), applying the method of separation of variables [43], The resulting differential operator C, in the 
transversal direction can be identified and applied to the corresponding basic equation [44] 

m = b(y) (27) 

and its Dirichlet boundary condition (BC) 

8(y m ) = gm- (28) 


The corresponding GF G(y,y) describes the solution for a differential equation with the same differential operator C in y and 
ErsatzBC the same boundary condition, but a special right-hand side inhomogeneity. This inhomogeneity is a Dirac impulse 
iny. The Dirac function is the neutral element of the convolution operation, thus the resulting solution of this problem can 
be used to find solutions for any other inhomogeneity. The basic form of the GF for a differential operator Cj of the form 


Cji 9 = 


1 _#_ 

aj dy 2 


,9. 


(29) 


reads 


G(y|y)= -aH(y-y)sinh(a(y-y))+C 1 (y)e“ 3 '4-C 2 (y)e ay , (30) 

with C] 2 depending on the boundary conditions g m . The index j denotes the equation the approach is applied to; in this case 
either v or K. Convoluting G(y|y) with the right-hand side of Eq. (27) in y and taking the BC into account yield 


■9(y) = G(y,y)*b(y) 

otj m °y 


(31) 


where the index m indicates a boundary. This solution is valid as long as aj is independent of y. This quantity may still vary 
along any other direction. 


4.1. Closure assumptions 

We restrict ourselves to the goal of finding an analytical solution for the problem considered. Due to the different nature 
of Eqs. (26b)-(26f) a coupled solution of the equations does not exist. However, solving the weakly coupled equations 
separately in segregated steps yields an insight as to which parts impede from finding an analytical solution. The definition 
of closure assumptions for such terms gives access to an analytical approach of the equations. Before the cross-sectional 
averaging procedure is carried out, the applied modeling is introduced. 
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The axial momentum equation (26c) can be rearranged such that the differential operator for the axial velocity iq 
becomes 



(32) 


Unfortunately, due to the convective term ugdu^/dx the operator £„[Ui] is not only an operator in they direction. Thus, an 
analytical GF cannot be found for this operator, only numerical solutions are available. This is in conflict with the aim of this 
paper to provide an analytical solution for the cross-sectional part of the set of governing (26). Therefore, the term indicated 
by A! either has to be modeled in terms of iq or formulated as inhomogeneity on the right-hand side with defined 
transversal shape (i.e. in this step p^). The latter approach implies that the term Aj is no longer a function of the axial 
acoustic velocity iq and thus does not influence the formulation of the GF any more. 

A similar convective term occurs in the fluid energy equation (26e). Substituting the acoustic density by Eq. (26a) 
facilitates writing the differential operator in terms of the fluctuating temperature T-i 



(33) 


This expression also includes a term (A 2 ), depending on the axial derivative of the parameter itself. Due to its similar 
nature, the modeling is chosen in accordance with A). As the transversal profile of the acoustic velocity iq is defined by the 
solution of Eq. (26c), all right-hand side terms of Eq. (26e) are known, and thus can be invoked in the analytical solution. 
Therefore, a reasonable modeling for both convective terms Ai and A 2 has to be found before the application of the GF 
method is carried out. Finding a model that represents the a priori unknown relevant physical impact of the convective term 
Ai is not straightforward. Iterative methods are often used to solve the problem of modeling an unknown term: Starting 
from a relative simple approach, the achieved solution can be used as a model in a following step. This procedure is repeated 
until the desired accuracy is reached. The convergence of this method strongly depends on the initial approach. Here, we 
present two alternative derivations. If higher accuracy is demanded, the solutions achieved may be used as a starting point 
for an iterative derivation. 

The first suggestion is a very general approach. Solving the equations without the bothersome terms A,- yields an 
analytical solution, which can be inserted in an iterative procedure. Especially the first iteration is advantageous for 
developing an understanding of the derivation process. Neglecting the convective terms reduces the length of the system of 
equations, which makes them manageable in a better way. Further, one may argue that at least for small convective scales in 
terms of Ma, the contribution of such components in an acoustically compact region should not affect the acoustic fields 
decisively. 

The second way is rather problem restrictive. As an initial solution, the zero mean flow solution is used. In contrast to a 
neglect of A it it describes a physical problem for special mean conditions. This idea leads to an immense increase in the 
complexity of the system. The number of emerging terms is many times that of the first approach. As this modeling is part of 
the first and the second of six derivation steps leading to the final set of ODEs, this increase causes a vast rise in the number 
of terms in the resulting set of equations. Thus, a simplified version of inserting the quiescent mean solution can be chosen 
by only considering the spatially averaged contribution of the terms <fyq /dx = (b(y/- l )/dx)\ Ua = 0 . This approach includes a more 
physical contribution of the term than completely neglecting it, like in the previous approach. However, this simplification is 
not motivated by physical considerations. One argument that may justify this step is the later cross-sectionally averaging of 
the entire system of equations, which makes it a look-ahead approach. The reader should be aware that the models 
presented here are intended to be used in the initial phase of an iterative derivation. Taking this averaging into account may 
lead to less accurate results of the first iteration step. Thus, this simplification only leads to slower convergence. 

In this publication, two of these modeling assumptions are carried out in order to observe their performance thoroughly: 
(i) In M I both convective terms are neglected, (ii) For M II the system of Eqs. (26a)-(26f) is solved for Ma = 0. This solution, 
where all gray boxed terms are zero, is inserted in a transversally averaged form. Table 3 summarizes the terms of both 
model assumptions. An iterative insertion of the corresponding solutions is beyond the scope of this publication. In the 
following sections, the modeled terms are replaced by their representative A,- indicating any contributions of the differently 
modeled source terms, which here are considered constant in the transversal direction. 

4.2. Cross-sectional averaging of the LNSEs 

Due to the large number of terms that occur in the spatial averaging process, the pertinent terminology is explained. The 
blackboard bold capital letters (e.g. E) indicate the equation, in which the replacement character is introduced at. The 
overbar (&) marks the terms stemming from transversally averaged equations. The indices refer to three different issues: (a) 
the corresponding prefactors (e.g. p^) the term multiplies with, (b) terms that are included in the integral (e.g. F„) or (c) the 
inclusion of a modeled term (A,). Table 4 lists the abbreviations and indices and gives reference to the corresponding 
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Table 3 

Overview over the modeled terms: For model M i, the convective terms are neglected. M II uses a spatially averaged, quiescent 
mean solution. 


Term 

Character 

Model 

M I 

M 11 

dll i 

A, 

0 

a<ui>| 

dX 



H 

ii 

o 

dT i 

A 2 

0 

^T,> 1 

dX 





Table 4 

Replacement characters and their indexing. 

Character 

Eq. 

Description 

Index 

Description 

E 

(26e) 

Energy conservation 

pi 

Acoustic pressure 

G 

(26a) 

Gas law 

pj 

Axial derivative of 

G 

(26a) 

Averaged gas law 

Ul 

Acoustic velocity 

Ik 

(26b) 

Intermediate averaged continuity 

u i 

Axial derivative of u A 




T, 

Acoustic temperature 

K 

(26b) 

Final averaged continuity 

Qi 

Acoustic density 

T 

(26e) 

Solution of energy 

Pi 

Axial derivative of q 1 


(26f) 

Conservations 

es 

Containing heat capacity ratio 

X 

(26c) 

Axial momentum 

F v 

Transversal viscous part 



conservation 

Fk 

Transversal thermal part 

A, 

(26c) 

Modeled momentum term 

A, 

Containing modeled term 1 

A2 

(26e) 

Modeled enthalpy term 

A2 

Containing modeled term 2 

equations in Appendix A.l. The averaging of the parameter Ui, Ti and <?] 

proceeds in the 

same order as in the zero mean 


flow derivation processed by in't panhuis [15], 

Like the mean transversal momentum Eq. (7d), the acoustic momentum conservation in y (26d) reveals that the acoustic 
pressure Pi(x,y) = p^(x) is constant over the cross-section of the channel. 

Applying the assumptions M I, M II defined in Section 4.1 into the momentum Eq. (26c) and using the differential 
operator of equation (29) allows one to rewrite it as 


— x Al + x Pl - 


dp i 
dX 


(34) 


u i I + i =0 


(35) 


where the coefficient a 2 in the differential operator C v [Ui] (see Eq. (31)) stands for 

du\ Wo 2 


a, =q \ IkJ - 

" ' dxj KjJ. 


(36) 


Arnott et al. [ 17] solved the GF problem of a similar differential operator, a constant inhomogeneity and a value of zero at the 
boundaries. Introducing the function derived in that work, 


F0) = 1 ~ 


cosh(gjy) 

cosh(gj) 


the solution of the transversal problem is obtained 

“1 = J C(y,y)<5(y-y) dy(x Al +x Pl ^ = F„(x Al +x Pl ^. 

Finally, we now perform the averaging in the y-direction. This yields 

<“i )= j o th dy= (l-/„)^X Al +X p ,-^. 
where the lowercase /denotes the Rott [12] function 

a i 

Eq. (39) is the first acoustic transport equation, which will be solved numerically later on. 


(37) 


(38) 


(39) 


(40) 
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Next, the density fluctuation q, is replaced in the energy equation (26e). Transforming it along with the solid energy 
equation (26f) into the same form as Eq. (34) and setting the coefficients that are constant along the transverse direction in 
the corresponding differential operator to 


. udT 
t,Cp l ~ lK+ T~dk 


2 _ i c s£?s 2Ta s 


2La / 

~i<K~ 


<p 2 Ks 

allows one to write them in combination with their coupling conditions as 


AdX l] — Ea 2 + E\ , ,f, F v + (E Pl + + (E P] - + E P| . f ;F^) 

Cs[Tt,sl = 0 
T'll + i = Ti,sl ± i =Tb 


dfh 

dX 


(41) 

(42) 

(43) 

(44) 

(45) 


I<- 


dy 


= K S 


dT i < 
dy 


(46) 


The solution of the differential equations with symmetry conditions at y=0 and y = 1+ F S /R leads to formulations of the 
temperature oscillations Ti, T l s that contain the coupling wall boundary temperature T b 


Tj = Ea 2 Fjc + E Ai ,f„Fk> + (E P] F k + E Pl>Fj ,Fj( j „)p 1 + (E Pi /FjH- E Pi ._f„Fk>)-|^--I-(1 -F K )T b 

Tv = (l-F K )T b . 


(47) 


(48) 


Inserting these equations into Eqs. (45) and (46) allows one to replace T l s and T b in the later formulation of the acoustic 
temperature transport. The constant terms in Ck[Ti] and dslT] s ] reveal again thermal Arnott functions, which only differ in 
the factor a.. The mixed term 


Fkv = 


-F K +yrF„ 


(49) 


-l + 'JSr 

in Eq. (47) is the result of the convolution of the GF and the viscous Arnott function. In this formulation, the square ratio of 
the operator constants a|/a^ was replaced by 'fir with respect to the implicit solution IMP. When mean flow is neglected, 
this ratio reduces to the Prandtl number Pr of the fluid. Next, the Neumann condition (46) can be applied to find 

//k> 


E A2 +E Al , f „ 


+ 


')BrE Pl + E Pl , Fi 


TlrE pp + E P] , a 


Ik 
’\ f K 

Ik,u 

Ik 


-1 

-1 

-1 


Pi 

dp i 

dX ’ 


whereas E es and e s are substitutes for 


E __ e 5_ 

W1+E S ) 


e S = i 


KWk 


(50) 


(51) 


(52) 


" KsPrfs 

Inserting Eq. (50) into Eq. (47) and introducing new variables T, E (see Appendix A.l) finally yields the expression for the 
temperature fluctuation 


Ti — T Al +Ta 2 f k F, ( 4- Ea, FK 4- ("E Pl + T Pl ,y K F k + E P| F kjj) p | 4- (T~ Pl ' T T Pl ' ,f k Fk 4- E P] ,F,F k,v) ^ . 


(53) 


This solution as well as its cross-sectional averaged form is now merged with the acoustic ideal gas law (26a) to obtain a 
formulation for the density fluctuation qi and its transversal average (<?,) 

dPi 


Pi = (Ta, +Ta 2 f k Fk + E Ai ^F k , 1 ,)G7' 1 + [G Pl + (T Pl 4-T Pl ,f k F k -|- E P]jP> ,F K _ i/ )G7 1 ]pi + (Tp 1 ' + T Pi ',f k Fk + E Pi ',f 1 .F)c,i*) | Gf 1 


ax 


(di) — G a ,a 2 T^PiFi + G P] 


dp i 
dX ’ 


(54) 

(55) 
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The average of the transversal acoustic component (Vj) vanishes in the integral form of the continuity (26b). Inserting the 
formulation for {q^) and d{Qr)/dx into the cross-sectional averaged conservation of mass (Eq. (26b)) yields 

= K Ai ,a 2 + Pi + ^PT~^r“*“ »W<Ui>. (56) 

Finally, using Eq. (39) leads to the second transport ODE 


d(Ul) 

dX 


A Pi' 




( 1 -/, 


<Ui> 


(57) 


that describes the axial transport of the acoustic velocity (iq) solely in terms of pressure fluctuation. The first two terms as 
well as the term X A , in the first transport ODE (Eq. (39)) show the impact of the closure terms. For M I the terms X A] and 
K A] , A2 vanish. If the closure assumptions A] 2 are formulated in terms of iq and pr (e.g. M II) these terms contribute to all 
components of the system matrix A of the transport equations (39) and (57) written as 


d^ 

dx 




(58) 


As only cross-sectional averaged acoustic quantities or quantities that do not depend on the transversal direction are 
considered from now on, the angular brackets are omitted. 


5. Spatial propagation and scattering matrices 

The accuracy of the closure assumptions chosen is validated against an existing one-dimensional solution denoted by 
IMP and a full CFD computation of one single pore. At first the x distributions for a single frequency are compared. Then, the 
frequency dependent scattering behavior is investigated in terms of scattering matrices. 


5.J. Scattering matrices 


In the investigated problem, the spatial distribution of p 1 (x) and iq(x) depends on the investigated frequency. Thus, a full 
characterization of the quality of closure terms by comparison of the spatial results for one frequency is impossible. 
Furthermore, the boundary conditions, i.e. the phase between the acoustic variables and their amplitude ratio strongly affect 
the solution. Hence, only one spatial profile is not enough to determine the full impact of mean flow. Therefore, a more 
general investigation tool has to be chosen: the acoustic scattering matrices. 

Standard low-order network models of wave propagation [45-47] based on one-dimensional acoustics use scattering 
matrices for the stability prediction of acoustic systems. Guedra et al. [48] used measured TA core data of this type to predict 
the onset of various TA engines. Scattering matrices describe the linear correlation of the characteristic wave amplitudes, 
which are often denoted as incoming and outgoing waves, in a causal sense. These complex valued amplitudes correlate to 
the acoustic variables by 

Pr=pc(f+G) (59) 


Ur =T-Q. 


(60) 


Here, the wave traveling in the positive axial direction is named T whereas the upstream traveling wave is designated Q. 

A variety of numerical and experimental identification methods of such matrices either apply the two-load-one-source 
location (see [49]) or the two-source location method introduced by Munjal et al. [50]. The common idea of both approaches 
is to create a set of two independent acoustic states. Inserting these into the definition of the scattering matrix 


hid tdd 
rim fdu 


(61) 


allows one to determine the linear coefficients for an acoustic element. Here, the first method is applied to the system by 
shifting the phase of the pressure boundary condition by ^/2. The integration of the ODEs is processed for both sets of BCs 
and a discrete number of frequencies. The same integration technique as for the investigation of the spatial distribution 
of pr and iq is applied in this case. 


5.2. Numerical integration of the thermoacoustic transport equations 

The frequency dependent boundary value problem is now solved numerically in the complex space using standard 
integration methods. Given the boundary conditions at the cold reference end, a standard numerical integration method can 
be used to solve the set of equations for the acoustic field (58) for a predefined mean temperature T(x ) and velocity 
distribution u(x). 

Here, a standard mathematics software tool is used to compute the relevant spatial distributions in segregated steps. 
At first, a shooting method 1 51 ] with an accuracy of 10~ 12 is used to solve Eq. (20) in combination with given mean 
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temperatures at both ends and a given cold side inlet Ma number. In a second step, quadratic interpolated values of the 
resolved temperature distribution T(x) are used to solve the ODE system for the acoustic field. This approximation is 
applicable, because the highest order spatial derivative of T occurring is d 2 T/dx 2 . An explicit linear multistep Adams- 
Bashforth method [51] of variable order and step size is applied. For quiescent conditions, the implementation was tested 
against the existing freeware DeltaEC [24,52] from the Los Alamos group around G.W. Swift. The code achieved deviations 
below 1 percent, which we ascribe to the different discretization and implementation of material properties. 


5.3. Validation cases 

Both closure approaches are compared for up to three different cases. One condition without mean flow Ma 0 = 0 and two 
conditions with mean flow at Mam 3 x 10~ 4 are investigated. For the latter two different Peclet numbers in a pore of 
e = 1 /52 are considered. First, a linear mean temperature profile is caused by a vanishing Pe number, which is denoted by 
Maii n . Then, convection is increased by setting Pem 30. This causes a clearly exponential-like distribution ( Ma exp ). Table 5 
summarizes the validation cases used. 

The presented system of thermoacoustic transport equations is validated against the existing implicit solution (IMP) of 
in't panhuis [37] and Swift [14] and full compressible CFD finite volume simulations. The same solution method is applied to 
the implicit and explicit one-dimensional ODEs for better comparability. The thermo-physical properties of the fluid are 
interpolated from tabulated values [53], 

At first, the spatial evolution at Wo = 1.74 is observed for Ma 0 and Ma exp using air at standard conditions. The changes in 
convective contribution to the mean heat flux from case Ma nn to Ma exp allow inferences on the importance of the 
temperature distribution in combination with mean flow. 

To compare the results to a full CFD simulation, the ID computation is adapted in two points. The boundary values for p, 
and u, are extracted from the periodic CFD result at x=0 and applied to the ID simulations. Further, the solid temperature is 
fixed, because the reference simulation cannot capture conjugate heat transfer. As a direct consequence of this restriction, 
the factor es is set to zero in the one-dimensional linear PDEs. 

The scattering matrices are computed for a Womersley number range of 0.02-5. The frequencies a> are determined such 
that they form an equally dense grid in the Wo space. For each frequency the system of equations (58) is solved once. From 
these data points the complex valued scattering matrix is constructed. 


5.3.1. Full CFD simulation 

The validation simulations are performed with the open source CFD code OpenFOAM [54], The full laminar compressible 
Navier-Stokes equations are solved for the domain sketched in Fig. 2. The geometrical scales are a length of L= 0.052 m and 
a half pore height of R = 1 x 10 -3 m. The wall temperatures are fixed in time. The simulation domain of the pore is extended 
by two free field areas with a length of 0.05 m, which are appended at both ends of the stack channel. The block structured 
grid consists of cells with a maximum axial size of Ax = 0.1 mm. Outside the pore, the grid is coarsened using a stretching 
factor of maximum 1.15 avoiding acoustic reflections at cell interfaces. The transversal direction is discretized by 10 cells. 
They experience an exponential refinement towards the solid wall boundary. Applying discretizing schemes that are of at 
least fourth order, this grid reproduces the smallest acoustic boundary layers considered here. Compared to the analytical 
solution [ 14] an overall accuracy of at least 99 percent was reached. Due to the fine meshing the acoustic Courant-Friedrich- 
Lewy (CFL) condition requires to use a time step of At= 1.25 x 10~ 7 s. 

The mean temperature at the wall is fixed such that the cross-sectionally averaged downstream temperature deviates up 
to 1 percent between the cases. The acoustic signals are excited at the upstream and downstream BCs located at both ends of 
the free fields. Those non-reflecting boundary conditions of the Euler characteristic BC [55] type use local one-dimensional 
inviscid (LODI) [56] relations. They are enhanced with a relaxation term that was proposed by Rudy and Strikwerda [57] and 
a plane wave masking filter technique [58], This combination showed a good non-reflecting behavior in different CFD codes 
[59,60], A time dependent term is further added for exciting the acoustic wave that is entering the domain. The amplitude is 
fixed such that the acoustic velocity u, reaches 5 percent of the characteristic velocity p/qc. 

The single frequency simulation needs to reach a quasi-steady state. Thus, it is processed until the acoustic variables 
Pi, ui.T, at three predefined plane locations (upstream end, downstream end and center of the pore) do not vary for more 
than 0.1 % 0 from the last three periods of oscillation. The time line is extracted with a sufficient temporal resolution 
(At = 200 /T). Averaging over five oscillation periods is performed for post-processing. 


Table 5 

Overview over the validation cases. 


Case 

Designation 

Mach number 

Peclet number 

Case 1 

Ma 0 

Ma = 0 

Pe = 0 

Case 2 

Ma,m 

Mam3x 10- 4 

Pe->0 

Case 3 

Mflexp 

Mfl«3x IQ- 4 

Pe & 30 
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Before the spatial averaging at the observation frequency is applied, the time domain information is transformed into 
frequency space and filtered. 

The multi-frequency simulations are performed applying pre-filtered random binary excitation signals. The CFD/SI 
method described in [61] is used for post-processing. This method and its successor the LES/S1 method (see [36,62,63]) 
apply the Wiener-filter method on S1SO (single input single output) and MISO (multiple input single output) systems to 
identify the frequency dependent transfer coefficients from time varying signals and responses. A laminar regime is 
observed, hence, no turbulent noise is induced into the simulation. The accuracy of such simulations scales with the signal 
length and the inverse of the frequency. The errors of results for wavelengths, fitting at least half inside the observed time 
line, are approximated to be in the range of + 1 percent. 


6. Results 


6.1. Spatial propagation for a single frequency 


The solid lines depicted in Fig. 4 represent the axial profiles of amplitude of the acoustic pressure (a) and the velocity 
fluctuation (b). Two cases are investigated: quiescent conditions (Ma 0 ) and the mean flow affected case Ma exp 3 x 10~ 4 . 
The displayed phase <P (dashed with empty symbols) is displayed relative to the reference phase of the acoustic pressure p^ 
located at the upstream end of the pore. 

All boundary conditions are chosen to be equal except for BC of the mean upstream velocity. Nevertheless, the upstream 
pressures (p, | x = 0 ) change with Mach number because of the different TA interaction processes and the imperfect 
reflectivity of the BCs. The variation of the axial evolution of |pj | for approximately 50 percent visualizes the influence of 
TA effects on the pressure field. Additionally, higher mean flow causes an increase in phase change of p ^ over the pore. The 
velocity distribution remains almost unaffected. This may be explained by the imposed characteristic waves and their phase 
difference of 0(p^)-0(u^)m> r/2 creating a standing wave-like condition near a velocity anti-node. 

Fig. 5 shows the relative magnitude and absolute phase error of the three investigated models (IMP, M I, M II) compared 
to the CFD solution for quiescent mean flow (Ma 0 = 0). The perfectly matching error curves indicate an identical behavior of 
all three tested quasi-one-dimensional models in comparison to the CFD simulation for the zero mean velocity condition. 
The maximum error in amplitude never exceeds 2.5 percent. The combination of this magnitude value and phase errors of 
up to 3° reveals a very good agreement of both simulation types. The differences are mainly ascribed to the different, multi¬ 
dimensional CFD conditions at the pore in- and outlet. The pressure amplitude and phase error exceeds the velocity devi¬ 
ation by a factor of two. 

In the presence of mean flow, the ID models do not behave identically (see Fig. 6). The ID acoustic pressure alternates 
around the CFD data. The decrease of the |Aiq| near the downstream position indicates a similar behavior at a larger length 
scale. The phase errors of all ID models behave equal: <P(Aiq) does not change with Mach number, while <£(Ap,) doubles. 
The amplitude errors of both transport variables p,, Ui also increase by this factor. The peak magnitude errors occur at 

0.8. Here, the explicit models allow a reduction of one fourth with respect to the IMP results. Model M II only slightly 
improves the results of M I. 

Altogether, Fig. 6 indicates that the newly derived models perform with higher accuracy than the implicit model IMP. 
This improvement especially concerns the axial distribution of the amplitude values, whereas the phase error of both iq, p, 
is almost unaffected by mean flow. Nevertheless, the TA mechanisms strongly depend on frequency. Thus, these single 
frequency observations do not imply a similar trend for all frequencies and mean flow conditions. Therefore, the prediction 
of the acoustic scattering behavior is investigated. 


(a) (b) 




Fig. 4. Axial profiles of amplitude (solid line, filled symbols) and phase (dashed line, empty symbols) of (a) acoustic pressure p i and (b) acoustic velocity Ui 
in a stack for two cases (Ma 0 in black and Ma exp with blue triangles) extracted from CFD. While the velocity Ui barely changes neither with space nor with 
mean flow, the impact of the latter is clearly visible for the pressure p lm (For interpretation of the references to color in this figure legend, the reader is 
referred to the web version of this paper.) 
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Fig. 5. Amplitude and phase error of the one-dimensional models IMP (•), MI(^) and M II (|) relative to the CFD results of Ma 0 (no mean flow, Pe=0). 
The filled symbols denote gain values and phase values are depicted with empty symbols. The acoustic pressure error is depicted in (a), (b) displays the 
velocity error. The errors from all three 1D cases behave identically. 
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Fig. 6. Error of the one-dimensional models IMP (•), MI (^) and MII (|) compared to the CFD results of Ma exp (mean flow with exponential temperature 
profile, Pe~30) for (a) acoustic pressure and (b) velocity. The filled symbols denote gain values and phase values are depicted with empty symbols. 
The amplitudes of the explicit models M I and M II deviate less from CFD than the implicit model IMP. 


6.2. Scattering matrices of the example pore 

First, the impact of mean flow and temperature field is investigated for the CFD reference data. Later, each case is 
compared to the results of the quasi-one-dimensional models IMP, M I and M II. 

Fig. 7 displays the CFD/SI results of the three cases selected. The indexing of the matrix elements denotes the 
contribution of the incoming (first index) to the outgoing (second index) characteristic wave amplitude. Furthermore, t and r 
indicate, whether the wave is transmitted or reflected respectively. Solid lines denote magnitude values of the complex 
valued coefficients scaling with the left y-axes, while the dashed lines mark the phases belonging to the right axes. 

The amplitudes of the reflection coefficients |r| in Fig. 7 are relatively unaffected by mean flow and smaller than those of 
the transmission coefficients |f|. Therefore, the small discrepancies in |r| are negligible and the variations in <£(r) play a 
minor role. 

The slope of the phase of the transmission coefficients 0{t) is constant for varying Mach and Peclet numbers. The wiggles 
at Wo > 2 are artifacts of the system identification procedure (see [63]). General considerations reveal that these oscillations 
are spurious. In this region the influence of the TA effect, which steadily decreases with frequency, becomes marginal. 

It was shown by Holzinger and Polifke [64] that the diagonal components of the scattering matrices control the TA 
energy transformation processes. Therefore, the amplitudes in the low Wo region are important for an accurate prediction of 
the scattering performance of a mean flow affected regenerator. The combination of a constant temperature gradient (blue 
curves in Fig. 7) and mean flow in the Ma lin case leaves the amplitude profile of t ud and t iu almost unaffected. There is only a 
slight reduction/enlargement visible for low frequencies. This corresponds to a decreasing conversion performance in such 
situations. 

In Fig. 8 the models IMP, M I and M II are validated against the CFD data for Ma 0 = 0. All three one-dimensional models 
overlap perfectly. This indicates an identical implementation for this case. The coefficients r uu , r dd match quite well for the 
CFD/SI data and the one-dimensional identification methods. On the contrary, the one-dimensional solution overestimates 
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Fig. 7. Scattering matrix of a slab pore at AT = 150 K and different Mach and Peclet numbers evaluated with the CFD/SI approach: Ma 0 (no mean flow, 
Pe=0), Ma lin (mean flow with linear temperature profile, Re . 0) and Ma exp (mean flow with exponential temperature profile. Re fw 30 ). Amplitude is 
depicted using a solid line and phase values are shown by dashed lines. The impact of the mean flow is strongest in both transmission coefficients of the 
scattering matrix at small Wo. 


the TA effect for low values of Wo, although the CFD data error underestimates its contribution in this range. The sample 
length of 1 s only allows an accurate identification up to Wo> 0.5. The CFD/SI method tends to predict the identity matrix for 
lower frequencies in transmission dominated acoustic elements. 

Fig. 9 compares the data for Ma\ in , i.e. Ma f» 3 x 10~ 4 and a vanishing Peclet number. Like the CFD data in Fig. 7, the 
explicit one-dimensional tools predict a slightly smaller TA interaction. Due to the identical mean temperature profile, the 
implicit model IMP is not at all affected by the velocity field and thus delivers results identical to Ma 0 conditions. Hence, the 
impact of including mean flow in the system of transport equations is small in the Ma Kn configuration. Thus, all models are 
applicable within a certain accuracy. 

The good agreement of the phases is also valid for all phase values of IMP and CFD data in the Pe f» 30 case shown in 
Fig. 10. Especially the phase of the reflected components <P(r) agree far better in the low frequency region. However, due to 
the low magnitude of these coefficients, the contribution to the overall scattering behavior is negligible. 

The phase differences of the ID results in the Wo < 1 region of t ud and t du are much smaller than the latter. The IMP 
formulation predicts an amplitude profile identical to Ma 0 for both transmission coefficients. They overestimate the 
influence of the TA effect by a value of 0.1, which is undesirable. Although the CFD prediction attenuates the amplification of 
the downstream (cold to hot) traveling wave by approximately 95 percent of the Ma 0 case, IMP does not consider this 
change. In contrast to IMP, M I predicts lower contributions of the thermo-viscous acoustic boundary layers to the 
transmission performance than the CFD/SI. When the frequency is reduced below Wo fa 1.5 the contributions strongly tend 
back to unity. In this range, the accordance of M II and the CFD simulation is very good. The curves of both identification 
processes almost match. Due to the steeper temperature gradient, the contributions of combined mean temperature and 
velocity terms mentioned for Pe-> 0 play a more profound role here. They capture the reduced energy transformation from 
heat and acoustics. It has to be stated that if Womersley numbers lower than 0.5 are considered, the solution becomes 
unphysical. The accumulation of numerical errors in the explicit mean velocity terms causes chaotic scattering predictions 
for the transmission coefficients that extend the displayed amplitude range by far [65], This problem is also observed for M I 
although occurring at far lower frequencies. 
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Fig. 8. Scattering matrix of a slab pore at AT = 150 K and Ma 0 (no mean flow, Pe=0). Amplitude is depicted using a solid line for the CFD or filled symbols 
in case of the one-dimensional models. Phase values are shown by dashed lines or empty symbols respectively. All one-dimensional models IMP (•), 
M I (|) predict equal values for the trivial case. They match quite well with the CFD solution. 


6.3. Physical explanations 

Finding a physical explanation for the thermoacoustic transport behavior under mean flow conditions is not 
straightforward. Two ideas are presented. While the former emphasizes the physics of the problem, the latter considers 
the mathematical aspects of the LNSEs. 

Reid [30] already proposed an explanation of the reduction in energy transmission by considering the temperature of gas 
parcels inside the stack. This temperature is displayed in Fig. 11 for stagnant and mean flow conditions respectively. Reid 
argues that the enclosed area (gray) in the T-x diagram scales with the energy transported by a fluid parcel. The area in the 
presence of mean flow should be smaller than for stagnant conditions. Thus, the thermoacoustic effect should be less 
effective with mean flow, which is indeed observed for the cases considered in this paper. At lower frequencies the ratio of 
the number of oscillations that a parcel experiences while passing through the stack to its time of residence is smaller, thus 
this effect would account for the reduction of the thermoacoustic hump for low Womersley numbers. 

The second, more mathematical approach is a consideration of the cross-sectional averaging functions occurring at the 
application of Green's function technique. The parameters aj occurring in the differential operator Cj of both the axial 
momentum equation (26c) and the energy equation (26e) are located on the positive imaginary axis of the complex plane. 
The additional terms in Eqs. (32) and (33) lead to a shift of aj away from this axis. The shape of the thermoacoustic hump in 
the scattering matrices correlates with the cross-sectional averaging function Fj of the pore considered [17,64]. If the 
argument of the function varies, its result also changes drastically. For low Wo, the Helmholtz number k is also small. Thus, 
the ratio of the new upcoming terms to the first term in the differential operator becomes smaller. Hence, the change in the 
resulting scattering matrix becomes larger for low Wo. This argument also explains the importance of the modeled 
convective terms Aj and A 2 , because they directly affect the cross-sectional shape function. They compensate the impact of 
the other terms such that the thermoacoustic hump is not reduced too much. 

The arguments presented provide some insight into the effects of mean flow on the scattering matrices of 
thermoacoustic stacks at high Peclet numbers Pe. However, it is admitted that the physical effects leading to the reduction 
of the thermoacoustic hump at low frequencies are not completely understood. In addition, we must concede that the 
second argument does not imply that the thermoacoustic hump should always be reduced in the presence of mean flow. 
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Fig. 9. Scattering matrix of a slab pore at AT = 150 K and Mai in (mean flow with linear temperature profile, Pe-»0). Amplitude is depicted using a solid line 
for the CFD or filled symbols in case of the one-dimensional models. Phase values are shown by dashed lines or empty symbols respectively. The one¬ 
dimensional models IMP (•}, M I f^j and M II (|) predict similar values for this case. The explicit models are slightly closer to the results from CFD/SI. 


These arguments only provide a small insight into the physics of the thermoacoustic effect in mean flow conditions. 
However, the one-dimensional modeling results match very well to the CFD data. This good agreement can be seen as an a 
posteriori justification for correctness of the presented model. Such methods are commonly applied when new modeling 
approaches are developed. 


7. Further improvements 

The system of equations (58) derived in Section 4 used the two simplest modeling approaches, i.e. neglecting the 
convective terms for M I and inserting a spatial averaged formulation for quiescent mean conditions for M II. Furthermore, 
the mean parameters were restricted to be constant in the transversal direction. This procedure kept the expressions short, 
but both simplifications may be replaced by more accurate methods. 

If laminar mean flow problems are considered, a parabolic velocity and temperature profile develops in the channel. 
Including such formulations into the LNSEs (26) leads not only into ay-dependency of the terms considered here, but also to 
the occurrence of additional terms that scale with transversal derivatives of the mean quantities. The first term in the 
differential operators C of Eqs. (32) and (33) then has a non-constant contribution. This differential operator leads to a 
combination of Weber and Gamma functions for the transversal momentum equation (26d). Additionally, the convective 
terms become functions of y and thus, the complexity of the system increases. Nevertheless, the derivation process still 
allows an analytical one-dimensional description of the acoustic transport process. 

The shift from transversally averaged model approaches for A, to its y-dependent equivalent may lead to more accurate results. 
As the GF of the problem does not change, the increase in computational costs is marginal. A third approach, which sequentially 
reinserts the solutions obtained by M I and M II for A,- also potentially improves the solutions, although the numerous terms 
inserted into the derivation at an early stage generate a huge number of expressions in the system matrix A(x). 

The ideas presented lead to a more physical description of the configuration and are subject to further investigations. 
Nevertheless, the further improvement in accuracy will probably not justify the increasing computational costs. This 
prediction is supported by the fact that the enhancement reached by M II increased the computational time by a factor 
of ten. 
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Fig. 10. Scattering matrix of a slab pore at AT = 150 K and Ma exp (mean flow with exponential temperature profile, Pe % 30). Amplitude is depicted using a 
solid line for the CFD or filled symbols in case of the one-dimensional models. Phase values are shown by dashed lines or empty symbols respectively. 
M I(^) and M II (|) equalize the prediction quality of IMP(») except for the phase of the reflection coefficients at Wo<1.5. In all other quantities the more 
complex model (M II) matches the CFD solution better than the implicit model (IMP). 



Fig. 11. Temperature of a fluid parcel considered by Reid [30]. The enclosed area (gray) of a mean flow affected parcel (b) is smaller than at stagnant mean 
conditions (a). 


8. Conclusions 

The one-dimensional explicit models derived in this study improve the accuracy of the description of the thermo-viscous 
interaction in narrow pores with mean flow. Earlier implicit models only incorporate the change in mean temperature 
distribution, which then interacts with the acoustic variables. The resulting scattering predictions of these implicit one¬ 
dimensional approaches are not affected by different, mean flow controlled temperature shapes. If mean flow is explicitly 
considered in the TA transport equations, the prediction of the scattering behavior is improved. An increase in complexity of 
the invoked closure assumptions leads to better prediction results for steep temperature profiles. Thus, the increased 
computational costs should be accepted in favor of a better predictability for Wo < 2. At very low frequencies the explicit 
models fail, because of the accumulation of numerical errors. 
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Although enhanced modeling approaches will be investigated in future, their practicability is questionable. This 
publication showed a good agreement of the reduced one-dimensional model M II to multi-dimensional CFD simulations 
in non-zero mean flow conditions. The expected gain in accuracy for such conditions occurring rarely in the literature is 
marginal. Time-averaged flow effects denoted as streaming, which are part of most TA devices, lead to non-parabolic mean 
flow profiles. If their magnitude is large, such conditions may necessitate a more detailed description of transversal 
components of the mean flow field, as well as transversally non-constant closure assumptions. If these enhancements are 
not sufficient, an iterative insertion of the obtained systems of equations may lead to satisfactory accuracy in the predicted 
results. All these improvements are accompanied by an increase in computational costs. At some point, they presumably do 
not justify the gain in accuracy any longer and the predictions may be carried out using coupled 1D-CFD simulations. 

Concluding, the usage of the enhanced quasi-one-dimensional model derived here allows accurate predictions for mean 
flow affected thermoacoustic devices. The development of such apparatus beneficially applying or at least handling inherent 
mean flow may lead to an increase in commercial application of the thermoacoustic effect. 

Appendix A. Abbreviations in formulas 

A. J. Axial momentum equation 


Comparing Eqs. (26c)—(34) under consideration of Eq. (32) directly leads to 

v 0~r)T 
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dx 
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~dX 
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A.2. Energy equation 

Inserting the solution of Eq. (34) into (26e) causes the abbreviations 
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Here, the terms E^ 2 and E Al vanish for 1V1 I. After combining both energy equations (26e) and (26f) the new occurring 
variables are 
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(A.3c) 

(A.3d) 

(A.3e) 

(A.3f) 


Again (A.3c) and (A.3f) vanish for the simpler model. 

A.3. Gas law 

The first insertion process of the resulting energy (53) and momentum equation (34) into the gas law (26a) yields two 
new parameters 

y 
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which are combined to 
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in the averaging process. Again the last term becomes zero for M I. 

A.4. Conservation of mass 

Finally, all derived formulas are inserted in the conservation of mass (26b) and the intermediate values 
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for the second transport PDE of the acoustic variables containing the term K Ai ,a 2 
model derivations. 
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